home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / tm_test.pro < prev    next >
Text File  |  1997-07-08  |  6KB  |  188 lines
  1. ;$Id: tm_test.pro,v 1.4 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       TM_TEST
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the Student's t-statistic and the probability
  11. ;       that two vectors of sampled data have significantly different means. 
  12. ;       The default assumption is that the data is drawn from populations with
  13. ;       the same true variance. This type of test is often refered to as the 
  14. ;       T-means Test.
  15. ;
  16. ; CATEGORY:
  17. ;       Statistics.
  18. ;
  19. ; CALLING SEQUENCE:
  20. ;       Result = TM_TEST(X, Y)
  21. ;
  22. ; INPUTS:
  23. ;       X:    An n-element vector of type integer, float or double.
  24. ;
  25. ;       Y:    An m-element vector of type integer, float or double.
  26. ;             If the PAIRED keyword is set, X and Y must have the same
  27. ;             number of elements.
  28. ;
  29. ; KEYWORD PARAMETERS:
  30. ;       PAIRED:   If set to a non-zero value, X and Y are assumed to be 
  31. ;                 paired samples and must have the same number of elements.
  32. ;
  33. ;       UNEQUAL:  If set to a non-zero value, X and Y are assumed to be from
  34. ;                 populations with unequal variances.
  35. ;
  36. ; EXAMPLE
  37. ;       Define two n-element vectors of tabulated data.
  38. ;         X = [257, 208, 296, 324, 240, 246, 267, 311, 324, 323, 263, 305, $
  39. ;               270, 260, 251, 275, 288, 242, 304, 267]
  40. ;         Y = [201, 56, 185, 221, 165, 161, 182, 239, 278, 243, 197, 271, $
  41. ;               214, 216, 175, 192, 208, 150, 281, 196]
  42. ;       Compute the Student's t-statistic and its significance assuming that 
  43. ;       X and Y belong to populations with the same true variance.
  44. ;       The result should be the two-element vector [5.5283890, 2.5245510e-06],
  45. ;       indicating that X and Y have significantly different means.
  46. ;         result = tm_test(X, Y)
  47. ;       
  48. ; PROCEDURE:
  49. ;       TM_TEST computes the t-statistic of X and Y as the ratio;
  50. ;       (difference of sample means) / (standard error of differences) and 
  51. ;       its significance (the probability that |t| could be at least as large
  52. ;       large as the computed statistic). X and Y may be of different lengths.
  53. ;       The result is a two-element vector containing the t-statistic and its
  54. ;       significance. The significance is a value in the interval [0.0, 1.0]; 
  55. ;       a small value (0.05 or 0.01) indicates that X and Y have significantly
  56. ;       different means. 
  57. ;
  58. ; REFERENCE:
  59. ;       Numerical Recipes, The Art of Scientific Computing (Second Edition)
  60. ;       Cambridge University Press
  61. ;       ISBN 0-521-43108-5
  62. ;
  63. ; MODIFICATION HISTORY:
  64. ;       Written by:  GGS, RSI, Aug 1994
  65. ;                    TM_TEST is based on the routines: ttest.c, tutest.c and
  66. ;                    tptest.c described in section 14.2 of Numerical Recipes, 
  67. ;                    The Art of Scientific Computing (Second Edition), and is 
  68. ;                    used by permission.
  69. ;-
  70.  
  71. function betacf, a, b, x
  72.   on_error, 2
  73.   eps   = 3.0e-7
  74.   fpmin = 1.0e-30
  75.   maxit = 100
  76.   qab = a + b
  77.   qap = a + 1.0
  78.   qam = a - 1.0
  79.     c = 1.0
  80.     d = 1.0 - qab * x / qap
  81.   if(abs(d) lt fpmin) then d = fpmin
  82.   d = 1.0 / d
  83.   h = d
  84.   for m = 1, maxit do begin
  85.     m2 = 2 * m
  86.     aa = m * (b - m) * x / ((qam + m2) * (a + m2))
  87.      d = 1.0 + aa*d
  88.      if(abs(d) lt fpmin) then d = fpmin
  89.      c = 1.0 + aa / c
  90.      if(abs(c) lt fpmin) then c = fpmin
  91.      d = 1.0 / d
  92.      h = h * d * c
  93.      aa = -(a + m) *(qab + m) * x/((a + m2) * (qap + m2))
  94.      d = 1.0 + aa * d
  95.      if(abs(d) lt fpmin) then d = fpmin
  96.      c = 1.0 + aa / c
  97.      if(abs(c) lt fpmin) then c = fpmin
  98.      d = 1.0 / d
  99.      del = d * c
  100.      h = h * del
  101.      if(abs(del - 1.0) lt eps) then return, h
  102.   endfor
  103.   message, 'Failed to converge within given parameters.'
  104. end
  105.  
  106. function gammln, xx
  107.   coff = [76.18009172947146d0,   -86.50532032941677d0,  $
  108.           24.01409824083091d0,    -1.231739572450155d0, $
  109.            0.1208650973866179d-2, -0.5395239384953d-5]
  110.   stp = 2.5066282746310005d0
  111.   x = xx
  112.   y = x
  113.   tmp = x + 5.5d0
  114.   tmp = (x + 0.5d0) * alog(tmp) - tmp
  115.   ser = 1.000000000190015d0
  116.   for j = 0, n_elements(coff)-1 do begin
  117.     y = y + 1.d0
  118.     ser = ser + coff[j] / y
  119.   endfor
  120.   return, tmp + alog(stp * ser / x)
  121. end
  122.  
  123. function ibeta, a, b, x
  124.   on_error, 2
  125.   if (x lt 0 or x gt 1) then message, $
  126.     'x must be in the interval [0, 1].'
  127.   if (x eq 0  or x eq 1) then bt = 0.0 $
  128.   else $
  129.     bt = exp(gammln(a + b) - gammln(a) - gammln(b) + $
  130.              a * alog(x) + b * alog(1.0 - x))
  131.   if(x lt (a + 1.0)/(a + b + 2.0)) then return, $
  132.     bt * betacf(a, b, x) / a $
  133.   else return, 1.0 - bt * betacf(b, a, 1.0-x) / b
  134. end
  135.  
  136. function tm_test, x0, x1, paired = paired, unequal = unequal
  137.  
  138.   on_error, 2
  139.  
  140.   if keyword_set(paired) ne 0 and keyword_set(unequal) ne 0 then $
  141.     message, 'Paired and Unequal keywords cannot be set simultaneously.'
  142.  
  143.   nx0 = n_elements(x0)
  144.   nx1 = n_elements(x1)
  145.  
  146.   if nx0 le 1 or nx1 le 1 then $
  147.     message, 'x0 and x1 must be vectors of length greater than one.'
  148.  
  149.   type = size(x0)
  150.  
  151.   if keyword_set(paired) ne 0 then begin
  152.     ;x0 and x1 are paired samples with corrected covariance.
  153.     if nx0 ne nx1 then message, $
  154.       'Paired keyword requires vectors of equal size.'
  155.     mv0 = moment(x0)
  156.     mv1 = moment(x1)
  157.     cov = total((x0 - mv0[0]) * (x1 - mv1[0]))
  158.     df = nx0 - 1
  159.     cov = cov / df
  160.     sd = sqrt((mv0[1] + mv1[1] - 2.0 * cov) / nx0)
  161.     t = (mv0[0] - mv1[0]) / sd
  162.     prob = ibeta(0.5*df, 0.5, df/(df+t^2))
  163.     if type[2] eq 4 then return, float([t, prob]) $
  164.     else return, [t, prob]
  165.   endif else if keyword_set(unequal) ne 0 then begin
  166.     ;x0 and x1 are assumed to have different population variances.
  167.     mv0 = moment(x0)
  168.     mv1 = moment(x1)
  169.     t = (mv0[0] - mv1[0]) / sqrt(mv0[1]/nx0 + mv1[1]/nx1)
  170.     df = (mv0[1]/nx0 + mv1[1]/nx1)^2 / $
  171.          ((mv0[1]/nx0)^2/(nx0 - 1.0) + (mv1[1]/nx1)^2/(nx1 - 1.0))
  172.     prob = ibeta(0.5*df, 0.5, df/(df+t^2))
  173.     if type[2] ne 5 then return, float([t, prob]) $
  174.     else return, [t, prob]
  175.   endif else begin
  176.     ;x0 and x1 are assumed to have the same population variance. 
  177.     mv0 = moment(x0)
  178.     mv1 = moment(x1)
  179.     df = nx0 + nx1 - 2
  180.     var = ((nx0 - 1)*mv0[1] + (nx1 - 1)*mv1[1]) / df 
  181.     t = (mv0[0] - mv1[0]) / sqrt(var*(1.0/nx0 + 1.0/nx1))
  182.     prob = ibeta(0.5*df, 0.5, df/(df+t^2))
  183.     if type[2] ne 5 then return, float([t, prob]) $
  184.     else return, [t, prob]
  185.   endelse
  186.  
  187. end
  188.